% $Id: parallel.tex 204 2009-04-20 01:59:36Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%

\mysection{Getting Started}
For the purposes of this section a cluster of computers is
simple a collection of computers joined by a network. They may be
something exotic like a Beowolf cluster with sophistcated
nertwork hardware or just a simple local network of machines
joined by one hub. As long as the amount of data that needs to be passed
from one computer to another is small one can successfully carry out
parallel computing even over a simple slow network.

At present the parallel distibuted computing module for \ADM\
uses the PVM (parallel virtual machine) software. 
To use the module you must install the PVM software on each computer 
in your cluster. There are installshield versions for MS operating systems
as well as various formats for Linux and other Unix systems.

Assuming that the software is installed on the various computers in your cluster
you must have an account on each system. ADM\ parallel processing
uses the master/slave model where one computer runs the master program and
there are one or more slave programs running on the computers in the
cluster. Note that more than one slave process can be run on each computer
and slave processes can be run on the same computer that is running the
master process.

\mysection{A simple example}
Assume that we have $n$ $m\times m$ matrices $M_i$ indexed by $i$.
Let the $m\times m$  matrix $S$ be given by
$$S=\sum_{i=1}^n (M_i+\lambda_i \hbox{\rm I})^{-1}$$
where I denoetes the $m\times m$ identity matrix.
The problem is to minimize the squared norm of $S$, that is the sum of the
squares of the elements of $S$ with respect to the $n$ parameters
$\lambda_i$.

The (unparallelized) code for the example is
\beginexample

DATA_SECTION
  int n
  int m
 !! m=10;
 !! n=100;
 !! random_number_generator rng(131);
  3darray C(1,n,1,m,1,m)
 !! for (int i=1;i<=n;i++) C(i).fill_randn(rng);
PARAMETER_SECTION
  init_vector lambda(1,m)
  init_number kludge(2)
  3darray M(1,n,1,m,1,m)
  3darray N(1,n,1,m,1,m)
  matrix S(1,m,1,m);
  objective_function_value f
PROCEDURE_SECTION
   int i;
   for (i=1;i<=n;i++)
   {
     M(i)=C(i);
     for (int j=1;j<=m;j++)
     {
       M(i)(j,j)+=lambda(i);
     }
   }
   // compute the inverse of each matrix. 
   // this is the code which will be parallelized
   for (i=1;i<=n;i++)
   {
     N(i)=inv(N(i));
   }

   S.initialize();
   for (i=1;i<=n;i++)
   {
     S+=N(i);
   }
   f=norm2(N);
   f+=.0001*norm2(lambda);
   f+=square(kludge);
  
TOP_OF_MAIN_SECTION
  
  long int arrmblsize=5000000;
\endexample
This example shows the simplest kind of code to parallelize.
There is an outer loop and the computations for each pass
through the loop are independent of each other.
\beginexample
   for (i=1;i<=n;i++)
   {
     N(i)=inv(N(i));
   }
\endexample
Here is the code for the master program for the parallelized computations.
\beginexample
DATA_SECTION
  int n
  int m
  ivector n_per_host(1,ad_nhost)
  ivector min_n(1,ad_nhost)
  ivector max_n(1,ad_nhost)
 !! m=10;
 !! n=100;
 !! random_number_generator rng(131);
  3darray C(1,n,1,m,1,m)
 !! for (int i=1;i<=n;i++) C(i).fill_randn(rng);
PARAMETER_SECTION
  init_vector x(1,m)
  init_number kludge(2)
  3darray M(1,n,1,m,1,m)
  3darray N(1,n,1,m,1,m)
  matrix S(1,m,1,m);
  objective_function_value f
PROCEDURE_SECTION
   int i;
   for (i=1;i<=n;i++)
   {
     M(i)=C(i);
     for (int j=1;j<=m;j++)
     {
       M(i)(j,j)+=x(i);
     }
   }
   determine_number_on_each_host();
   send_data(); // send M to slaves
   receive_data(); // get N back
   
   S.initialize();
   for (i=1;i<=n;i++)
   {
     S+=N(i);
   }
   f=norm2(N);
   f+=.0001*norm2(x);
   f+=square(kludge);
  
FUNCTION  determine_number_on_each_host
  int i;
  n_per_host=n/ad_nhost;
  int remainder=n%ad_nhost;
  for (i=1;i<=remainder;i++) n_per_host(i)++;


  min_n(1)=1;
  max_n(1)=n_per_host(1);
  for (i=2;i<=ad_nhost;i++) 
  {
    min_n(i)=min_n(i-1)+n_per_host(i-1);
    max_n(i)=min_n(i)+n_per_host(i)-1;
  }

FUNCTION send_data
  int i;
  for (i=ad_nhost;i>=1;i--) 
  {
    int bufid = pvm_initsend( PvmDataRaw );
    pvm_pack(gradient_structure::no_derivatives); // tell slaves if derivatives are calculated
    pvm_pack(M.sub(min_n(i),max_n(i)));
    adpvm_send((*ad_stid)[i-1],1);
  }

FUNCTION receive_data
  int i;
  cout << "starting pvm_receive" << endl;
  
  for (i=ad_nhost;i>=1;i--) 
  {
    cout << "starting pvm_receivei for host "<< i << endl;
    adpvm_recv((*ad_stid)[i-1],2);  // get the values from slave
    pvm_unpack(N.sub(min_n(i),max_n(i)));
    adpvm_end_receive();  // have got all the values
  }
 
GLOBALS_SECTION
  extern "C" {
  #include "pvm3.h"
  }
  #include <admodel.h>
  #include <adpvm.h>

  #if defined(WIN32)
  #  include <windows.h>
  #endif

  struct pvmhostinfo *ad_hostp=NULL;
  int ad_nhost=0;
  ivector * ad_stid=NULL;

TOP_OF_MAIN_SECTION
  
  long int arrmblsize=5000000;
  cout << "calling load library" << endl;
  if (load_adpvm_library() < 0)
  {
    cerr << "error loading pvm library" << endl;
    exit(1);
  }
  int result, check, i, narch;
  pvm_setopt(PvmRoute, PvmRouteDirect);  /* channel for communication */
  /* get and display configuration of the parallel machine */
  pvm_config( &ad_nhost, &narch, &ad_hostp );  /* get configuration */
  printf("I found the following hosts in your virtual machine\n");
  for (i = 0; i < ad_nhost; i++)
  { 
     printf("\\t%s\n", ad_hostp[i].hi_name);
  }
  ad_stid= new ivector(0,ad_nhost);
  
  for (i=0; i<ad_nhost; i++)				/* spawn processes on */			
  {						/* all physical machines */
    check=pvm_spawn("pslave", 0,PvmTaskHost,ad_hostp[i].hi_name, 1, &((*ad_stid)(i)));
    if (!check) 
       printf("Couldn't start process on %s\n",ad_hostp[i].hi_name); 
    else
       printf("Started process on %s\n",ad_hostp[i].hi_name); 
  }

\endexample
The code for sending data to the slave processes and getting it back has been 
put into two functions.aHere is the code to send the data. At this point
{\tt ad\_nhost} is the number of machines in the cluster and there is one slave
process running on each member of the cluster.

\beginexample
FUNCTION send_data
  int i;
  for (i=ad_nhost;i>=1;i--) 
  {
    int bufid = pvm_initsend( PvmDataRaw );
    pvm_pack(gradient_structure::no_derivatives); // tell slaves if derivatives are calculated
    pvm_pack(M.sub(min_n(i),max_n(i)));
    adpvm_send((*ad_stid)[i-1],1);
  }
\endexample
To begin the process of sending data one first calls the  {\tt pvm\_initsend}
function. This initializes a data buffer to hold the data prior
to sending them. The {\tt PvmDataRaw} option means that the data are not tranlated
into a format reconizeable my computers that use incompatible data
formats but are sent ``as is''. This is the faster option.
The line {\tt  pvm\_pack(gradient\_structure::no\_derivatives);}
should always be included first. This is because the slave program
has no way of knowing whether derivative calculations have been turned
on or off. It uses the value of this variable to control
derivative calculations.
The line {\tt pvm\_pack(M.sub(min\_n(i),max\_n(i)));} put the {\tt sub dvar3\_array}
which is determined by the integers {\tt min\_(i)} and {\tt max\_(i)} into the buffer.
the {\tt sub} member function has been defined to produce sub objects of the
\ADM\ container objects. This is very convenient for communication in
parallelized code.
Now we are ready to send the data. This is accomplished with the
    {\tt adpvm\_send((*ad\_stid)[i-1],1);} 
  function. Here {\tt ad\_stid} is a
pointer to an {\tt ivector} which contains the identifier for
the i'th slave process in its i-1'th component.

The code for the slave process is not an \ADM\ program. It is a C++ progam which
uses the Autodif libraries which underly \ADM\ as well as the ADPVM module.
\beginexample

#include <stdio.h>
extern "C" {
#include <pvm3.h>
}
#include <admodel.h>
#include <adpvm.h>
#if !defined(linux)
#  include <windows.h>
#endif
main(int argc,char * argv[])
{
  gradient_structure gs(25000000);
  time_t now;
  char name[12], buf[60];
  int ptid;

  load_adpvm_library();
  ppvm_setopt(PvmRoute, PvmRouteDirect);
  
  do
  {
    dvar3_array N;

    // ***************  begin receive block ***********************************
     int ptid=ppvm_parent();
     adpvm_recv(ptid,-1);
     pvm_unpack(gradient_structure::no_derivatives);  // see if derivatives are
     pvm_unpack(N);  // see if derivatives are
     adpvm_end_receive();
    // ***************  end receive block ***********************************
  
     int nmin=N.indexmin();
     int nmax=N.indexmax();

     for (int i=nmin;i<=nmax;i++)
     {
       N(i)=inv(N(i));
     }
  
    // ***************  begin send block ***********************************
    ppvm_initsend(PvmDataDefault);/* allocate message buffer */	
    pvm_pack(N);
    adpvm_send(ptid, 2);/* send buffer to master */
    // ***************  end send block ***********************************

    slave_gradcalc();
  }
  while(1);
}

\endexample

The code to receive the data from the parent program is
\beginexample
     int ptid=ppvm_parent();
     adpvm_recv(ptid,-1);
     pvm_unpack(gradient_structure::no_derivatives);  // see if derivatives are
     pvm_unpack(N);  // see if derivatives are
     adpvm_end_receive();
\endexample
The first line {\tt int ptid=ppvm\_parent();} get a handle to the parent process
so that the slave process knows where to get its data from. The line
     { \tt adpvm\_recv(ptid,-1);} receives the data and puts it into a buffer.
The line
    {\tt  pvm\_unpack(gradient\_structure::no\_derivatives);} put the value of the 
the flag into {\tt  gradient\_structure::no\_derivatives);} which turns
derivative calculations off or on as desired. 
The line {\tt  pvm\_unpack(N);}
puts the values of the {\tt dvar3\_array} into the {\tt dvar3\_array N}.
If {\tt N} is unallocated this function will allocate the right sizes for 
accepting the data that are sent. This means that one does not need to
know or calculate in advance the shape of the container object. One can simple
declare a {\tt dvar3\_array N} and everything works.
Finally every receive block which includes variable objects must be
terminated with the line {\tt adpvm\_end\_receive();}. If you leave this out
the code will not work.

\mysection{Improving the performance of the  code}
The above code was run with n=100 and m=20. The unparallelized code
was much faster when using a simple 10MBit/sec hub. However it
is possible to speed up the code by being a bit smarter.
Since {\tt C} is constant we only need to pass it to the
slave processes once.the size of {\tt C} is 100$\times$20$\times$20$\times$8 bytes
or 320,000 bytes. Also once the inverses of the matrices have been calculated
by the slave  processes the inverse matrices can be added so that only
the sum need to be passed back to the master process.
Here is the code for parallel1.tpl which is modified to take advantage
of these optimzations. 
\beginexample
DATA_SECTION
  int n
  int m
  ivector n_per_host(1,ad_nhost)
  ivector min_n(1,ad_nhost)
  ivector max_n(1,ad_nhost)
 !! m=20;
 !! n=100;
 !! random_number_generator rng(131);
  3darray C(1,n,1,m,1,m)
 LOC_CALCS
  for (int i=1;i<=n;i++) C(i).fill_randn(rng);
  int i;
  n_per_host=n/ad_nhost;
  int remainder=n%ad_nhost;
  for (i=1;i<=remainder;i++) n_per_host(i)++;
  min_n(1)=1;
  max_n(1)=n_per_host(1);
  for (i=2;i<=ad_nhost;i++) 
  {
    min_n(i)=min_n(i-1)+n_per_host(i-1);
    max_n(i)=min_n(i)+n_per_host(i)-1;
  }
  ad_constant_flag=1;
  for (i=ad_nhost;i>=1;i--) 
  {
    int bufid = pvm_initsend( PvmDataRaw );
    pvm_pack(m); // send C to slaves
    pvm_pack(C.sub(min_n(i),max_n(i))); // send C to slaves
    adpvm_send((*ad_stid)[i-1],1);
  }
  ad_constant_flag=0;
 END_CALCS
PARAMETER_SECTION
  init_vector x(1,n)
  vector f0(1,ad_nhost)
  3darray M(1,n,1,m,1,m)
  3darray sumN(1,ad_nhost,1,m,1,m)
  matrix S(1,m,1,m);
  objective_function_value f
PROCEDURE_SECTION
   int i;
   send_data(); // send x to slaves
   receive_data(); // get sumN back
   
   S.initialize();
   for (i=1;i<=ad_nhost;i++)
   {
     S+=sumN(i);
   }
   f=norm2(S);
   f+=.001*norm2(x);
  
FUNCTION send_data
  int i;
  for (i=ad_nhost;i>=1;i--) 
  {
    int bufid = pvm_initsend( PvmDataRaw );
    pvm_pack(gradient_structure::no_derivatives); // tell slaves if derivatives are calculated
    pvm_pack(x(min_n(i),max_n(i))); // send x to slaves
    adpvm_send((*ad_stid)[i-1],1);
  }

FUNCTION receive_data
  int i;
  
  for (i=ad_nhost;i>=1;i--) 
  {
    adpvm_recv((*ad_stid)[i-1],2);  // get the values from slave
    pvm_unpack(sumN(i));
    adpvm_end_receive();  // have got all the values
  }
 
GLOBALS_SECTION

  extern int ad_constant_flag;

  extern "C" {
  #include "pvm3.h"
  }
  #include <admodel.h>
  #include <adpvm.h>

  #if defined(WIN32)
  #  include <windows.h>
  #endif

  struct pvmhostinfo *ad_hostp=NULL;
  int ad_nhost=0;
  ivector * ad_stid=NULL;

TOP_OF_MAIN_SECTION
  
  long int arrmblsize=5000000;
  cout << "calling load library" << endl;
  if (load_adpvm_library() < 0)
  {
    cerr << "error loading pvm library" << endl;
    exit(1);
  }
  int result, check, i, narch;
  pvm_setopt(PvmRoute, PvmRouteDirect);  /* channel for communication */
  /* get and display configuration of the parallel machine */
  pvm_config( &ad_nhost, &narch, &ad_hostp );  /* get configuration */
  printf("I found the following hosts in your virtual machine\n");
  for (i = 0; i < ad_nhost; i++)
  { 
     printf("\\t%s\\n", ad_hostp[i].hi_name);
  }
  ad_stid= new ivector(0,ad_nhost);
  
  for (i=0; i<ad_nhost; i++)				/* spawn processes on */			
  {						/* all physical machines */
    check=pvm_spawn("pslave1", 0,PvmTaskHost /* | PvmTaskDebug */ ,ad_hostp[i].hi_name, 1, &((*ad_stid)(i)));
    if (!check) 
       printf("Couldn't start process on %s\\n",ad_hostp[i].hi_name); 
    else
       printf("Started process on %s\\n",ad_hostp[i].hi_name); 
  }

\endexample
The process of sending data to the slaves has been broken up into two parts. Once the
{\tt d3\_array C} has been calculated a send block is set up and the parts of
{\tt C} are shipped off to the slave processes. At present the program can not
detect that derivatives are not being calculated here.  To generate the
correct code the line {\tt ad\_constant\_flag=1;} must be inserted into the code and
{\tt extern int ad\_constant\_flag;} must be declared in the globals section.
At the end of the send block the line {\tt ad\_constant\_flag=0;}  turns on the
generation of the proper derivative code.
\beginexample
  ad_constant_flag=1;
  for (i=ad_nhost;i>=1;i--) 
  {
    int bufid = pvm_initsend( PvmDataRaw );
    pvm_pack(m); // send m to slaves
    pvm_pack(C.sub(min_n(i),max_n(i))); // send C to slaves
    adpvm_send((*ad_stid)[i-1],1);
  }
  ad_constant_flag=0;
\endexample
The code for the slave processes is
\beginexample
#include <stdio.h>
extern "C" {
#include <pvm3.h>
}
#include <admodel.h>
#include <adpvm.h>
#if !defined(linux)
#  include <windows.h>
#endif
#define  __declspec(dllexport) 

extern int ad_constant_flag;

main(int argc,char * argv[])
{
  gradient_structure gs(25000000);
  time_t now;
  char name[12], buf[60];
  int ptid;

  load_adpvm_library();
  pvm_setopt(PvmRoute, PvmRouteDirect);
  
  d3_array C;
  int m;
  // ***************  begin receive block ***********************************
   ad_constant_flag=1;
   ptid=pvm_parent();
   adpvm_recv(ptid,-1);
   pvm_unpack(m);  // see if derivatives are
   cout << "m = " << m << endl;
   pvm_unpack(C); 
     adpvm_end_receive();
   ad_constant_flag=0;
  // ***************  end receive block ***********************************
  do
  {
    dvar_vector x;

    // ***************  begin receive block ***********************************
   ptid=pvm_parent();
     adpvm_recv(ptid,-1);
     pvm_unpack(gradient_structure::no_derivatives);  // see if derivatives are
     pvm_unpack(x);  // see if derivatives are
     adpvm_end_receive();
    // ***************  end receive block ***********************************
  
     int nmin=C.indexmin();
     int nmax=C.indexmax();
     dvar_matrix M(1,m,1,m);
     dvar_matrix N(1,m,1,m);
     N.initialize();

     dvariable f0=0.0;
     for (int i=nmin;i<=nmax;i++)
     {
       M=C(i);
       for (int j=1;j<=m;j++)
       {
         M(j,j)+=x(i);
       }
       N+=inv(M);
     }
  
    // ***************  begin send block ***********************************
   ptid=pvm_parent();
    cout << " starting send " << endl;
    pvm_initsend(PvmDataDefault);/* allocate message buffer */	
    pvm_pack(N);
    adpvm_send(ptid, 2);/* send buffer to master */
    cout << " finished send " << endl;
    // ***************  end send block ***********************************

    cout << "starting slave gradcalc " << endl;
    slave_gradcalc();
    cout << "finished slave gradcalc " << endl;
  }
  while(1);
}

\endexample
The code for the extra receive block is below. As in the send block for the
master process.  The {\tt ad\_constant\_flag} must be set so that the proper code is
produced.
\beginexample
   ad_constant_flag=1;
   ptid=pvm_parent();
   adpvm_recv(ptid,-1);
   pvm_unpack(m); 
   pvm_unpack(C); 
   adpvm_end_receive();
   ad_constant_flag=0;
\endexample
With these modifications the parallelized code ran about twice as fast
the unparallelized code.
